/*
**    Copyright (C) 2003-2011 Institute for Systems Biology
**                            Seattle, Washington, USA.
**
**    This library is free software; you can redistribute it and/or
**    modify it under the terms of the GNU Lesser General Public
**    License as published by the Free Software Foundation; either
**    version 3 of the License, or (at your option) any later version.
**
**    This library is distributed in the hope that it will be useful,
**    but WITHOUT ANY WARRANTY; without even the implied warranty of
**    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
**    Lesser General Public License for more details.
**
**    You should have received a copy of the GNU Lesser General Public
**    License along with this library; if not, write to the Free Software
**    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
*/

package org.systemsbiology.stats;

import org.apache.commons.math.random.MersenneTwister;
import org.apache.commons.math.random.RandomData;
import org.apache.commons.math.random.RandomDataImpl;

import java.io.*;
import java.util.Arrays;

/**
 * Created by IntelliJ IDEA.
 * User: anorberg
 * Date: 10/18/11
 * Time: 1:57 PM
 *
 * Class to generate, save, and load a lookup table of result values
 */
public class HartiganDipLut implements Serializable{
    static final long serialVersionUID = 943L;

    public static final int DEFAULT_TRIALS = 1000001; //one million and one
    public static final int MANY_TRIALS = 10000001; //ten million and one

    private static int[] dup(int[] in){
        int[] ret = new int[in.length];
        System.arraycopy(in, 0, ret, 0, in.length);
        return ret;
    }

    private static double[] dup(double[] in){
        double[] ret = new double[in.length];
        System.arraycopy(in, 0, ret, 0, in.length);
        return ret;
    }

    private static final int[] SUPERDENSE_N_SAMPLES = {
            4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
            20, 25, 30, 35, 40, 45, 50,
            60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200,
            225, 250, 275, 300, 325, 350, 375, 400, 425, 450, 475, 500,
            550, 600, 650, 700, 750, 800, 850, 900, 950, 1000,
            1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000,
            2250, 2500, 2750, 3000, 3250, 3500, 3750, 4000, 4250, 4500, 4750, 5000,
            5500, 6000, 6500, 7000, 7500, 8000, 8500, 9000, 9500, 10000
    };
    private static final int[] DEFAULT_N_SAMPLES = {
            4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
            20, 25, 30,
            40, 50,
            75, 100,
            150, 200, 250, 300, 350, 400, 450, 500,
            600, 700, 800, 900, 1000,
            1250, 1500, 1750, 2000,
            3500, 4000, 4500, 5000
    };
    private static final double[] DEFAULT_P_VALUES = {
            0,
            0.01, 0.02, 0.05,
            0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
            0.95, 0.98, 0.99,
            0.995, 0.998, 0.999,
            0.9995, 0.9998, 0.9999,
            0.99995, 0.99998, 0.99999,
            1
    };
    private static final double[] HIGH_PRECISION_P_VALUES = {
            0,
            0.00001, 0.00002, 0.00005,
            0.0001, 0.0002, 0.0005,
            0.001, 0.002, 0.005,
            0.01, 0.02, 0.05,
            0.1, 0.15, 0.2, 0.25, 0.3, 0.35,
            0.4, 0.5, 0.6,
            0.7, 0.75, 0.8, 0.85, 0.9,
            0.95, 0.98, 0.99,
            0.995, 0.998, 0.999,
            0.9995, 0.9998, 0.9999,
            0.99995, 0.99998, 0.99999,
            0.999995, 0.999998, 0.999999,
            1
    };

    public int[] getSuperdenseNSamples(){
        return dup(SUPERDENSE_N_SAMPLES);
    }
    public int[] getDefaultNSamples(){
        return dup(DEFAULT_N_SAMPLES);
    }
    public double[] getDefaultPValues(){
        return dup(DEFAULT_P_VALUES);
    }
    public double[] getHighPrecisionPValues(){
        return dup(HIGH_PRECISION_P_VALUES);
    }

    private void fill(double[] n, RandomData rand){
        for(int k = 0; k < n.length; ++k){
            n[k] = rand.nextUniform(0, 1);
        }
    }

    private double[] p_values;
    private int[] n_samples;
    private double[][] calculated_values; //index 0: n_samples. index 1: p_values.

    /**
     * Print a debug message if and only if debugging is enabled. Convenience method that takes
     * a whole bunch of hairy-looking if statements out of the real code. If the "loud" flag is not
     * set, this method is a no-op.
     *
     * @param loud Whether or not to print anything at all.
     * @param message Message to print, if and only if the loud flag was set.
     */
    private static void spew(boolean loud, String message){
        if(loud){
            System.out.println(message);
        }
    }

    public HartiganDipLut(int[] nIn, double[] pIn, int precision, boolean loud){
        spew(loud, "Initializing LUT...");
        p_values = new double[pIn.length];
        System.arraycopy(pIn, 0, p_values, 0, pIn.length);
        n_samples = new int[nIn.length];
        System.arraycopy(nIn, 0, n_samples, 0, nIn.length);

        calculated_values = new double[n_samples.length][];
        RandomData rand = new RandomDataImpl(new MersenneTwister());

        int percentIncr = precision / 100;
        int[] dexes = new int[p_values.length];
        for(int k = 0; k < dexes.length; ++k){
            dexes[k] = (int)((precision - 1) * p_values[k]);
        }

        spew(loud, "Beginning calculation.");
        double[] scratch = new double[precision];
        for(int nDex = 0; nDex < n_samples.length; ++nDex){
            int n = n_samples[nDex];
            spew(loud, "Begin size " + n);
            double in[] = new double[n];
            for(int trial = 0; trial < precision; ++trial){
                fill(in, rand);
                scratch[trial] = HartiganDipTest.dip(in).getDip();
                if(trial % percentIncr == 0){
                    spew(loud, "Size " + n + " now " + (100.0 * trial / precision) + "% complete.");
                }
            }
            spew(loud, "Sorting " + precision + " results for " + n + ".");
            Arrays.sort(scratch);
            double[] result = new double[p_values.length];
            for(int pDex = 0; pDex < p_values.length; ++pDex){
                result[pDex] = scratch[dexes[pDex]];
            }
            calculated_values[nDex] = result;
            spew(loud, "Row " + nDex + " calculated.");
        }
        spew(loud, "LUT built.");

    }

    public HartiganDipLut(int[] n, double[] p, int precision){
        this(n, p, precision, false);
    }

    public double interpolateUnimodalP(double dip, int n){
        if(n <= 3){
            return 1.0; //all tiny samples are purely unimodal
        }
        int nDex = Arrays.binarySearch(n_samples, n);
        double plusRatio = 0;
        int n0 = -1;
        int n1 = -1;
        int nDexNext = nDex + 1;
        if(nDex < 0){
            //miss
            nDexNext = -nDex - 1;
            nDex = nDexNext - 1;
        } else {
            //direct hit
            nDexNext = nDex + 1; //semi-moot
        }
        //check for being off the chart either way, then find N position
        if(nDex < 0){
            //n underflow. implies our table doesn't go low enough, since it has to cover 4 and values under 3
            //already got a 1
            throw new AssertionError("Lookup table not configured for small number of points: " + n + ". Always start the table at 4");
        } else if(nDexNext >= n_samples.length){
            //past or at the end
            //overlimit case: assume final row is asymptotic
            nDex = nDexNext = n_samples.length -1; //last index cram case
            n0 = n1 = n_samples[nDex];
            //plusRatio remains 0
        } else {
            //typical case: mid-table. Possibly a direct hit, doesn't matter at this point
            n0 = n_samples[nDex];
            n1 = n_samples[nDexNext];
            double miss = n - n0; //type upcast for math reasons
            double range = n1 - n0;
            plusRatio = miss / range;
        }

        double[] interpolatedLut = new double[p_values.length];

        double sqrtN0 = Math.sqrt(n0);
        double sqrtN1 = Math.sqrt(n1);

        for(int k = 0; k < p_values.length; ++k){
            double y0 = sqrtN0 * calculated_values[nDex][k];
            double y1 = sqrtN1 * calculated_values[nDexNext][k];
            interpolatedLut[k] = y0 + ((y1 - y0) * plusRatio);
        }

        //interpolatedLut is expected to also be sorted as a linearly weighted sum of two sorted lists.

        double target = Math.sqrt(n) * dip;

        int targetDex = Arrays.binarySearch(interpolatedLut, target);
        if(targetDex >= 0){
            //direct hit, but how broad is it?
            int low = targetDex;
            int high = targetDex;
            while(low >= 1){
                if(interpolatedLut[low-1] == target){
                    --low;
                } else {
                    break;
                }
            }

            while(high < interpolatedLut.length - 1){
                if(interpolatedLut[high+1] == target){
                    ++high;
                } else {
                    break;
                }
            }
            return 1.0 - ((p_values[low] + p_values[high])/2);
        } else {
            int nextOver = -targetDex - 1;
            if(nextOver <= 0){
                return 1.0; //off the end
            }
            if(nextOver >= p_values.length){
                return 0.0; //off the other end
            }
            //normal case: somewhere between
            double offset = target - interpolatedLut[nextOver - 1];
            double range = interpolatedLut[nextOver] - interpolatedLut[nextOver - 1];
            double ratio = offset/range;

            return 1.0 - (
                    p_values[nextOver - 1] + (
                            ratio * (p_values[nextOver] - p_values[nextOver - 1])
                            )
                    );
        }
    }

    public double calculateDipUnimodalP(double[] x){
        return interpolateUnimodalP(HartiganDipTest.dip(x).getDip(), x.length);
    }

    public static void main(String[] args) throws IOException {
        if(args.length != 3){
            System.err.println("Wrong number of arguments.");
            System.err.println("Required: [precise|fast] precision outfile");
            System.exit(1);
        }
        boolean precise = args[0].equalsIgnoreCase("precise");
        if(!precise){
            if(!args[0].equalsIgnoreCase("fast")){
                System.err.println("Unknown mode: " + args[0]);
                System.err.println("Legal modes: precise, fast");
                System.exit(2);
            }
        }
        int precision = -1;
        try{
            precision = Integer.parseInt(args[1]);
        } catch(NumberFormatException fail){
            System.err.println(args[1] + " is not an identifiable int.");
        }
        if(precision <= 0){
            System.out.println("Using default precision");
        }

        ObjectOutputStream objWrite = new ObjectOutputStream(new FileOutputStream(args[2]));
        HartiganDipLut calculatedLut = null;
        System.out.println("Beginning calculation...");
        if(precise){
            calculatedLut = new HartiganDipLut(
                    HartiganDipLut.SUPERDENSE_N_SAMPLES,
                    HartiganDipLut.HIGH_PRECISION_P_VALUES,
                    precision <= 0 ? HartiganDipLut.MANY_TRIALS : precision,
                    true) ;
        } else {
            calculatedLut = new HartiganDipLut(
                    HartiganDipLut.DEFAULT_N_SAMPLES,
                    HartiganDipLut.DEFAULT_P_VALUES,
                    precision <= 0 ? HartiganDipLut.DEFAULT_TRIALS : precision,
                    true);
        }
        System.out.println("Calculated.");
        objWrite.writeObject(calculatedLut);
        System.out.println("Saved.");
    }
}
